
    #Article:    Legislature Resizing with Rent-Seeking Politicians: The Impact of Executive-Legislature Coalitions			
    #Author:     Anderson Frey <anderson.frey@rochester.edu>			    

    # This script produces all Tables and Figures in the Online Appendix
    # Tables are produced in a way that they can be directly copied/pasted in a LaTex document
    
    # It needs the following files:  data_main.csv  | data_reelection.csv 
    
    # Please upload all libraries and functions below, and run the setup before attempting to produce any Tables or Figures

    #Set the working directory to the folder that contains the replication files, use the function below
    
    #setwd()


    # Libraries ========================================================================================================
    
    
      library(data.table) 
      library(dplyr) 
      library(lfe)
      library(rdrobust)
      library(ggplot2) 
      library(rddensity)
      library(gridExtra)
      library(fastDummies)
      library(tidyr)
      library(boot)
      library(ggpubr)
      
  
    #===================================================================================================================

    
    # Functions ========================================================================================================
    
    
    # RD Plots
    .RDPlot <- function(data,cut,band,bins,degree,x.title,y.title,title,con,sz,brk,limi1,limi2) {
      
      if(degree == 1){mm <- (var ~ Score2 )}
      if(degree == 2){mm <- (var ~ Score2 + I(Score2**2) )}
      if(degree == 3){mm <- (var ~ Score2 + I(Score2**2) + I(Score2**3) )}
      
      data <- subset(data, !is.na(var))
      
      data$Score2 <- data$run - cut 
      data <- subset(data, abs(Score2) <= band)
      vec <- seq(cut, band, length.out =  (bins+1))
      vec_mid <- vec[-length(vec)] + diff(vec)/2
      data$count <- 1 ;     i <- 0
      
      
      for (iii in 0:1){
        cc <- iii 
        if(iii == 0) cc <- -1
        aux <- subset(data, t == iii)
        aux$groups <- as.factor(cut(cc*aux$Score2,breaks=vec, include.lowest = TRUE))
        R <- group_by(aux,groups) %>% summarise (var = mean(var), count = sum(count) ) %>% complete(groups)
        #R <- subset(R, !is.na(Var))
        R$t <- iii
        R$Score2 <- cc*vec_mid
        
        mdl <- lm(mm , data = aux) 
        pred <- predict(mdl, R, interval = "confidence", level = con, na.action = na.pass) 
        PRE <- pred[,1] ; R$pre <- PRE
        predSE <- predict(mdl, R, interval = "confidence", level =con, na.action = na.pass, se=TRUE)$se.fit
        t <- ((pred[1,1]-pred[1,2])/predSE[1])
        V <- as.matrix(vcov(mdl) )  
        
        if(degree == 1){X <- data.frame(1,R$Score2)}
        if(degree == 2){X <- data.frame(1,R$Score2,(R$Score2^2))}
        if(degree == 3){X <- data.frame(1,R$Score2,(R$Score2^2),(R$Score2^3))}
        
        X <- as.matrix(X) ; se_robust <- sqrt(diag(X %*% V %*% t(X)))
        LOW <- pred[,1] - t*se_robust ; HIGH <- pred[,1] + t*se_robust
        R$low <- LOW ;    R$high <- HIGH 
        
        if(iii == 0){R0 <- R ; plot_data <- R}else{R1 <- R ;plot_data <- rbind(plot_data,R)}
      }
      
      plot_data <- data.frame(plot_data)
      plot_data <- filter(plot_data, !is.na(count))
      plot_data$count2 <- plot_data$count/sum(plot_data$count)
      if(sz==0)  plot_data$count2 <- 1
      ymx <- max(plot_data$var)
      ymn <- min(plot_data$var)
      
      ggplot(data=plot_data, aes(x=Score2, y=var, ymax = ymx, ymin = ymn))   +   theme_bw() +
        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
        ylim(limi1,limi2)+
        geom_point( data=R0, aes(x=Score2, y=var, size=count), shape=21,  color = "gray25", bg="gray45") + #aes(size = count2/4)
        geom_point( data=R1, aes(x=Score2, y=var,size=count), shape=21,   color = "gray55", bg="gray75") +
        geom_vline(xintercept = 0, colour="gray25", linetype = "longdash") +
        geom_line(data=R1, aes(x=Score2, y=pre),colour="gray55",  linewidth=2, alpha=.95) +
        geom_line(data=R0, aes(x=Score2, y=pre),colour="gray25",  linewidth=2, alpha=.95) +
        ggtitle(title) + 
        theme(plot.title = element_text(family = "A", size=12))+
        xlab(x.title) + ylab(y.title) +
        theme(axis.text = element_text(colour="black", size=12, family="A")) +
        theme(axis.title = element_text(colour="black", size=12, family="A")) +
        geom_ribbon(data=R0, aes(ymin=low,ymax=high),alpha=0.10) +
        geom_ribbon(data=R1, aes(ymin=low,ymax=high),alpha=0.10) +
        scale_x_continuous(breaks=brk) +
        theme(legend.position="none")
      
      
    }
    
    # Regression output in LaTex format
    .LatexTable <- function(data,rownames,spaces){
      cols <- ncol(data)*2+1
      res <- data.frame(matrix(0,nrow(data),cols))
      res[,] <- "&"
      fill <- seq(from=2,to=cols,by=2)
      res[,fill] <- data
      res <- cbind(rownames,res)
      lastone <- rep("\tabularnewline",nrow(res))
      lastone[spaces] <- "\tabularnewline[0.2cm]"
      lastone[nrow(res)] <- "\tabularnewline[0.2cm]"
      res[,ncol(res)] <- lastone
      names(res) <- NULL
      return(res)
    }
    
    # Plot of marginal effects
    .NewPlot <- function(data,r,v,xvar,minD,maxD,bins1,bins2,shift,bb,x.title,y.title,
                         p.title,degree,multi,oneside,bb3, pt){
      
      
      Di <- matrix(0,(bins1),3)
      Di[,1] <- seq(from=minD, to=maxD, length.out=(bins1))
      label <- paste("t",xvar,sep=":")
      
      for (i in 1:nrow(Di)){
        Di[i,2] <- r["t",1] + r[label,1]*Di[i,1] 
        Di[i,3] <-(  v["t","t"] + (Di[i,1]**2)*v[label,label] + 2*Di[i,1]*v["t",label] )^.5
      }
      
      MT <- data.frame(Di) ; names(MT) <- c("D","coeff","sd")
      MT$h <- MT$coeff+1.96*MT$sd
      MT$l <- MT$coeff-1.96*MT$sd  
      
      jazz <- data.frame(matrix(0,(bins2),2))
      xx <- seq(from=min(MT$D), to=(max(MT$D)), length.out=(bins2+1))
      
      data <- data.frame(data)
      data$x <- data[,match(xvar,names(data))]
      
      i <- 2
      
      for (i in 1:nrow(jazz)){
        top <- xx[i+1]
        bot <- xx[i]
        if(is.na(top))  top <- xx[i+1] + 1
        jj <- subset(data, x >= bot & x < top )
        if(oneside==1) jj <- subset(data, x >= bot & x < top & t == 1)
        if(oneside==1) {datax <- subset(data,  t == 1)}else{datax <- data}
        jazz[i,1] <- (xx[i]+xx[i+1])/2
        
        jazz[i,2] <- nrow(jj)/nrow(datax)
      }
      names(jazz) <- c("D","dens")
      
      bb2 <- bb+shift
      
      if(pt==1){
        ggplot(data=MT, aes(x=D, y=coeff, ymin=0))   +    theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_line(data=MT, aes(x=D, y=coeff+shift),colour="black",  linewidth=1, alpha=.75) +
          theme(plot.title = element_text(family = "A", size=13))+
          geom_hline(yintercept=(0+shift), lty=3)+
          xlab(x.title) + ylab(y.title) +
          theme(axis.text = element_text(colour="black", size=13, family="A")) +
          theme(axis.title = element_text(colour="black", size=13, family="A")) +
          geom_line(data=MT, aes(x=D, y=h+shift),colour="black",  linetype="dashed",linewidth=.75, alpha=.75) +
          geom_line(data=MT, aes(x=D, y=l+shift),colour="black",  linetype="dashed",linewidth=0.75, alpha=.75) +
          ggtitle(p.title) +
          geom_bar(data=jazz, aes(x=D, y=dens*multi, group = 1), stat="identity", position="dodge",
                   size =0.1,color="gray35",fill="gray85", alpha=0.3) +
          scale_y_continuous(breaks=bb2, labels=bb )+
          scale_x_continuous(breaks=bb3)+
          coord_cartesian(ylim=c(bb2[1],bb2[length(bb2)]))
      }else{
        ggplot(data=MT, aes(x=D, y=coeff, ymin=0))   +    theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_line(data=MT, aes(x=D, y=coeff+shift),colour="black",  linewidth=1, alpha=.75) +
          theme(plot.title = element_text(family = "A", size=13))+
          geom_hline(yintercept=(0+shift), lty=3)+
          xlab(x.title) + ylab(y.title) +
          theme(axis.text = element_text(colour="black", size=13, family="A")) +
          theme(axis.title = element_text(colour="black", size=13, family="A")) +
          geom_line(data=MT, aes(x=D, y=h+shift),colour="black",  linetype="dashed",linewidth=.75, alpha=.75) +
          geom_line(data=MT, aes(x=D, y=l+shift),colour="black",  linetype="dashed",linewidth=0.75, alpha=.75) +
          geom_bar(data=jazz, aes(x=D, y=dens*multi, group = 1), stat="identity", position="dodge",
                   size =0.1,color="gray35",fill="gray85", alpha=0.3) +
          scale_y_continuous(breaks=bb2, labels=bb )+
          scale_x_continuous(breaks=bb3)+
          coord_cartesian(ylim=c(bb2[1],bb2[length(bb2)]))
      }
      
    }
    
    
    # For bootstrapping difference in coefficients from 2 RDD regressions
    .DiffRDD <- function( base, indices, b1,b2, controls, variable){
      
      library(fastDummies)
      library(rdrobust)
      library(dplyr)
      library(data.table)
      
      bb <- data.table( base[indices,] )
      ddd <- bb %>% select((fe))
      datafex <- dummy_cols(ddd)
      datafex <- datafex[,3:ncol(datafex)]
      
      
      bb$i <- bb %>% select(all_of(variable))
      
      bb$ti <- bb$t*bb$i
      bb$tir <- bb$t*bb$i*bb$run
      bb$ir <- bb$i*bb$run
      dataconx <- bb %>% select(all_of(controls))
      cv <- cbind(bb$i,bb$ti,bb$tir,bb$ir,dataconx,datafex)
      
      reg <- rdrobust(bb$var,bb$run,  b=b2, h=b1,  covs=cv,  masspoints="off") 
      a <- reg$coef[1]
      bb$i <- 1-bb$i
      
      bb$ti <- bb$t*bb$i
      bb$tir <- bb$t*bb$i*bb$run
      bb$ir <- bb$i*bb$run
      
      cv <- cbind(bb$i,bb$ti,bb$tir,bb$ir,dataconx,datafex)
      
      reg <- rdrobust(bb$var,bb$run,  b=b2, h=b1, covs=cv, masspoints="off") 
      b <- reg$coef[1]
      
      return(b-a)
      
    }
    
    
    #===================================================================================================================
    
    
    # Setup ============================================================================================================
    
      # set fonts for plots
      windowsFonts( A=windowsFont("Liberation Sans")) # Font for ggplot2 
      
      
      # Open file
      base <- read.csv("data_main.csv")
      
      
      # Adjust covariates
      base$lpoverty <- log(base$poverty)
      base$lvoters <- log(base$present08)
      base$lArea <- log(base$Area)
      base$rent.inc <- 0
      base$rent.inc[base$m.party.last %in% c(15,11,22,14)] <- 1
      base$pt_base <- ifelse(base$m.party %in% c(11,12,13,14,15,22,40),1,0)
      base$v <- (base$vote_win*base$last.w + base$op.vote_win*base$last.l )*100 / base$tv_winner 
      base$run <- base$residual*100 / base$tv_winner
      base$run <- (2*base$last.w-1)*base$run
      
      # Define treatment dummy and fixed-effects
      base$t <- base$last.w
      base$fe <- base$seats08
      
      # Parameters for optimal bandwidth algorithm 
      btype <- "cercomb1"
      qq <- 5
      vart <- "hc3"
      sdd <- sd(abs(base$run))
      
      
      # List of covariates
      controls <- c("qe","lgdppc","lbpc","lvoters","lpoverty","rent.inc","lArea","pt_base")
    
    
    
    #===================================================================================================================
    

    # Table A.1 ========================================================================================================
    
        # keeps only observations where the close election was: coalition vs. opposition
        new <- filter(base, MainSample==1)    

        ct <- 9 
        new$cap <- ct
        cuts <- c(15,30,50,80,120,160,300,450,600,750,900,1050,1200,1350,1500,1800,2400,3000,4000,5000,6000,7000,8000)
        for (i in 1:length(cuts)){
          new$cap[new$poplaw >= cuts[i]] <- ct+2
          ct <- ct+2
        }
        
   
        size12 <- c(0,15,30,50,80,120,160,300,450,600,750,900,1050,1200,1350,1500,1800,2400,3000,4000,5000,6000,7000,8000,90000000)
        
        lm12 <- data.table(matrix(0,length(size12)-1,3))
        
        count <- 9
        
        for (i in 1:nrow(lm12)){
          
          lm12[i,1] <- size12[i]
          lm12[i,2] <- size12[i+1]
          lm12[i,3] <- count 
          count <- count+2
        }
        
        size08 <- c(0,47619,95238,142857,190476,238095,285714,333333,380952,428571,476190,523809,571429,1000000,1121905,1243903,1365854,1487805,
                    1609756,1731707,1853658,1975609,4999999,5119047,5238094,5357141,5476188,5595235,5714282,5833329,5952376,6071423,
                    6190470,6309517,6428564,6547611,100000000)/1000
        
        lm08 <- data.table(matrix(0,length(size08)-1,3))
        
        count <- 9
        
        for (i in 1:nrow(lm08)){
          
          if(i==14) count <- 33
          lm08[i,1] <- size08[i]
          lm08[i,2] <- size08[i+1]
          lm08[i,3] <- count 
          count <- count+1
        }
        
        
        jazz <- data.table(c(size08,size12))
        jazz <- jazz %>% arrange(V1)
        jazz <- jazz[2:(nrow(jazz)-1),]
        
        names(lm12) <- c("low","high","cap")
        names(lm08) <- c("low","high","seats")
        
        names(jazz) <- "low"
        
        jazz <- left_join(jazz, lm12, by="low")
        jazz$high <- NULL
        for (i in 1:10){
          for (j in 1:nrow(jazz)){
            if(is.na(jazz$cap[j] )) jazz$cap[j] <- jazz$cap[j-1] 
          }
          
        }
        
        jazz <- left_join(jazz, lm08,by="low")
        jazz$high <- NULL
        for (i in 1:10){
          for (j in 1:nrow(jazz)){
            if(is.na(jazz$seats[j] )) jazz$seats[j] <- jazz$seats[j-1] 
          }
          
        }
        
        jazz$high <- 0
        jazz$high[1:(nrow(jazz)-1)] <- jazz$low[2:(nrow(jazz))]
        
        jazz <- jazz[1:(nrow(jazz)-1),]
        jazz$expansion <- jazz$cap-jazz$seats
        jazz$expansion2 <- round( jazz$expansion / jazz$seats , 3)*100
        
        jazz$high <- round(jazz$high,1)
        jazz$low <- round(jazz$low,1)
        tab <- jazz %>% select(low,high,seats,cap,expansion,expansion2)
        
        spaces <- nrow(tab)
        nam <- rep(" ",nrow(tab))
        jazz <- .LatexTable(tab,nam,spaces)
        
        
        # First Page
        print(jazz[1:30,3:14], row.names=FALSE)
        
        jazz[59,5] <- "-"
        # Second Page
        print(jazz[31:59,3:14], row.names=FALSE)
        
  
        rm(tab,jazz,lm08,lm12,size12,size08,ct,new,cuts,count)
    
    #===================================================================================================================
    

    # Table A.2 ========================================================================================================
   
        
        # keeps only observations where the close election was: coalition vs. opposition
        new <- data.table(filter(base, MainSample==1)   ) 
        
     
        
        # list of covariate names
        controls_n <- c("Votes per seat ('000)","GDP pc","Municipal Budget pc","Voters","Poverty",
                        "Municipal Area (km2)",
                        "Rent-Seeking Right (Mayor)", "PT Federal Base (Mayor)")
        
        

        
        for (out in 1:length(controls)){
          
          new$var <- new %>% select(controls[out])
          dataother <- new %>% select(controls[-out])
          res <- data.frame(matrix(0,3,5))
          
          ddd <- new %>% select(fe)
          datafe <- dummy_cols(ddd)
          datafe <- datafe[,3:ncol(datafe)]
          datacon <- new %>% select(all_of(controls))
          ddd <- new %>% select(reg)
          datasta <- dummy_cols(ddd)
          datasta <- datasta[,3:ncol(datasta)]
          
          cv <- datafe
          b <- rdbwselect(new$var, new$run,  p=1, q=qq, bwselect = btype)$bws ; b
          b1 <- b[1]
          b2 <- b1*(b[3]/b[1])
          
          for (i in 1:3){
            
            cv <- datafe
            if(i == 2) cv <- cbind(datafe,dataother)
            if(i == 3) cv <- cbind(datafe,dataother,datasta)
            
            reg <- rdrobust(new$var,new$run, covs=cv, b=b2, h=b1, q=qq, vce=vart)
            reg2 <- rdrobust(new$var,new$run, covs=cv, b=b2, h=b1, q=qq, vce=vart, level=90 )
            
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],2)
            ci2 <- round(reg$ci[3,2],2)
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            
            res[1,i+1] <- paste(format(co,nsmall=3),stars,sep="")
            res[2,i+1] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
            res[3,i+1] <- paste( "[" ,format(ci1,nsmall=2),",",format(ci2,nsmall=2), "]",sep="")
            if(i==3) res[1,i+2] <- format(round(b1/sdd,2),nsmall=2)
            
            
          }
          res[,1] <- controls[out]
          if(out==1){resx <- res}else{resx <- rbind(resx,res)}
          print(out)
        }
        
        resx <- resx[,2:5]
        resx[resx==0] <- ""
        
        spaces <- seq(from=3, by=3, length=length(controls))
        for (i in 1:length(controls_n)){
          nam1 <- c(controls_n[i]," "," ")
          if(i == 1){nam <- nam1}else{nam <- c(nam,nam1)}
        }
        tt <- .LatexTable(resx,nam,spaces)
        print(tt, row.names=F)
        
        
    #===================================================================================================================
        
        
    # Table A.3 ========================================================================================================
        
        # keeps only observations where the close election was: coalition vs. opposition
        data.table(filter(base, MainSample==1)   )   
        
        major <- 0
       
        if(major==0) ccc <- new %>% select(gender.0,college.0,new.0,incum.0,right.0,ptbase.0)
        if(major==1) ccc <- new %>% select(gender.1,college.1,new.1,incum.1,right.1,ptbase.1)
        namesccc <- names(ccc)
        
        # runs the regressions
        for (out in 1:ncol(ccc)){
          
          new$var <- new %>% select(namesccc[out])
          res <- data.frame(matrix(0,3,4))
          
          ddd <- new %>% select(fe)
          datafe <- dummy_cols(ddd)
          datafe <- datafe[,3:ncol(datafe)]
          datacon <- new %>% select(all_of(controls))
          ddd <- new %>% select(reg)
          datasta <- dummy_cols(ddd)
          datasta <- datasta[,3:ncol(datasta)]
          
          b <- rdbwselect(new$var, new$run, p=1, q=qq, bwselect = btype)$bws ; b
          
          b1 <- b[1]
          b2 <- b1*(b[3]/b[1])
          
          
          for (i in 1:3){
            
            
            cv <- datafe
            if(i == 2) cv <- cbind(datafe,datacon)
            if(i == 3) cv <- cbind(datafe,datacon,datasta)
            
            reg <- rdrobust(new$var, new$run,  covs=cv, b=b2, h=b1, q=qq, vce=vart) #; summary(reg)
            reg2 <- rdrobust(new$var,new$run, covs=cv, b=b2, h=b1, q=qq, level=90, vce=vart)
            
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],2)
            ci2 <- round(reg$ci[3,2],2)
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            
            res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
            res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
            res[3,i] <- paste( "[" ,format(ci1,nsmall=2),",",format(ci2,nsmall=2), "]",sep="")
            if(i==3) res[1,i+2] <- format(round(b1/sdd,2),nsmall=2)
            
            
          }
          
          if(out==1){resx <- res}else{resx <- rbind(resx,res)}
          print(out)
        }
        
        resx[is.na( resx)] <- ""
        resx$X4 <- NULL
        
        # produces the LaTex Table
        controls2_n <- c("Gender","Education","Newcomer","Incumbent","Rent-seeking Right", "PT Federal Base")
        
        
        spaces2 <- seq(from=3, by=3, length=length(controls2_n))
        for (i in 1:length(controls2_n)){
          nam1 <- c(controls2_n[i]," "," ")
          if(i == 1){nam <- nam1}else{nam <- c(nam,nam1)}
        }
        tty <- .LatexTable(resx,nam,spaces2)
        print(tty, row.names=F)
        
        
    #===================================================================================================================
        
        
    # Table A.4 ========================================================================================================
        
        # keeps only observations where the close election was: coalition vs. opposition
        data.table(filter(base, MainSample==1)   ) 
        
        major <- 1
        
        if(major==0) ccc <- new %>% select(gender.0,college.0,new.0,incum.0,right.0,ptbase.0)
        if(major==1) ccc <- new %>% select(gender.1,college.1,new.1,incum.1,right.1,ptbase.1)
        namesccc <- names(ccc)
        
        # runs the regressions
        for (out in 1:ncol(ccc)){
          
          new$var <- new %>% select(namesccc[out])
          res <- data.frame(matrix(0,3,4))
          
          ddd <- new %>% select(fe)
          datafe <- dummy_cols(ddd)
          datafe <- datafe[,3:ncol(datafe)]
          datacon <- new %>% select(all_of(controls))
          ddd <- new %>% select(reg)
          datasta <- dummy_cols(ddd)
          datasta <- datasta[,3:ncol(datasta)]
          
          b <- rdbwselect(new$var, new$run, p=1, q=qq, bwselect = btype)$bws ; b
          
          b1 <- b[1]
          b2 <- b1*(b[3]/b[1])
          
          
          for (i in 1:3){
            
            
            cv <- datafe
            if(i == 2) cv <- cbind(datafe,datacon)
            if(i == 3) cv <- cbind(datafe,datacon,datasta)
            
            reg <- rdrobust(new$var, new$run,  covs=cv, b=b2, h=b1, q=qq , vce=vart) #; summary(reg)
            reg2 <- rdrobust(new$var,new$run, covs=cv, b=b2, h=b1, q=qq, level=90, vce=vart)
            
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],2)
            ci2 <- round(reg$ci[3,2],2)
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            
            res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
            res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
            res[3,i] <- paste( "[" ,format(ci1,nsmall=2),",",format(ci2,nsmall=2), "]",sep="")
            if(i==3) res[1,i+2] <- format(round(b1/sdd,2),nsmall=2)
            
            
          }
          
          if(out==1){resx <- res}else{resx <- rbind(resx,res)}
          print(out)
        }
        
        resx[is.na( resx)] <- ""
        resx$X4 <- NULL
        
        # produces the LaTex Table
        controls2_n <- c("Gender","Education","Newcomer","Incumbent","Rent-seeking Right", "PT Federal Base")
        
        
        spaces2 <- seq(from=3, by=3, length=length(controls2_n))
        for (i in 1:length(controls2_n)){
          nam1 <- c(controls2_n[i]," "," ")
          if(i == 1){nam <- nam1}else{nam <- c(nam,nam1)}
        }
        tty <- .LatexTable(resx,nam,spaces2)
        print(tty, row.names=F)
        
        
        
    #===================================================================================================================
        
        
    # Table A.5 ========================================================================================================
        
        # keeps only observations where the close election was: coalition vs. opposition
        new <- data.table( filter(base, MainSample==1 ) )   
        
        
        new$lp <- as.numeric(substring(new$no_win,1,2))*new$last.w + 
          as.numeric(substring(new$op.no_win,1,2))*new$last.l
        new$rent.council <- ifelse(new$lp %in% c(15,11,22,14),1,0)
        
        
        interact <- c("patr","share","rent.council","v") 
        interact_n <- c("Patronage (binary)","Patronage (continuous)","Rent-seeking Right","Councilor's Safety") 
        
        ### Keep only data with donor information
        new$lp <- as.numeric(substring(new$no_win,1,2))*new$last.w + 
          as.numeric(substring(new$op.no_win,1,2))*new$last.l
        
        new$rent.council <- ifelse(new$lp %in% c(15,11,22,14),1,0)
        
        out <- 1
        for (out in 1:4){
          
          if (out <= 2){
            newx <- filter(new, !is.na(share) ) 
            newx$patr <- ifelse(newx$share > median(newx$share),1,0)
          }else{
            newx <- new
          }
          
          newx$var <- newx %>% select(interact[out])
          res <- data.frame(matrix(0,3,4))
          
          ddd <- newx %>% select(fe)
          datafex <- dummy_cols(ddd)
          datafex <- datafex[,3:ncol(datafe)]
          dataconx <- newx %>% select(all_of(controls))
          ddd <- newx %>% select(reg)
          datastax <- dummy_cols(ddd)
          datastax <- datastax[,3:ncol(datasta)]
          
          b <- rdbwselect(newx$var, newx$run, p=1, q=qq, bwselect = btype)$bws ; b
           
          b1 <- b[1]
          b2 <- b1*(b[3]/b[1])
          
          
          for (i in 1:3){
            
            cv <- datafex
            if(i == 2) cv <- cbind(datafex,dataconx)
            if(i == 3) cv <- cbind(datafex,dataconx,datastax)
            
            b <- rdbwselect(newx$var, newx$run, p=1, q=qq, bwselect = btype)$bws ; b
            
            b1 <- b[1]
            b2 <- b1*(b[3]/b[1])
            
            reg <- rdrobust(newx$var,newx$run, covs=cv, b=b2, h=b1, q=qq, vce=vart)
            reg2 <- rdrobust(newx$var,newx$run, covs=cv, b=b2, h=b1, q=qq, level=90, vce=vart)
            
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],2)
            ci2 <- round(reg$ci[3,2],2)
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            
            res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
            res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
            res[3,i] <- paste( "[" ,format(ci1,nsmall=2),",",format(ci2,nsmall=2), "]",sep="")
            if(i==3) res[1,i+1] <- format(round(b1/sdd,2),nsmall=2)
            
            
          }
          
          if(out==1){resx <- res}else{resx <- rbind(resx,res)}
          print(out)
        }
        resx[resx==0] <- ""
        
        spaces <- c(3,6,9,12)
        for (i in 1:length(interact)){
          nam1 <- c(interact_n[i]," "," ")
          if(i == 1){nam <- nam1}else{nam <- c(nam,nam1)}
        }
        tt <- .LatexTable(resx,nam,spaces)
        print(tt, row.names=F)
        
    #===================================================================================================================
        
        
    # Table A.6 ========================================================================================================
        
        a <- base
        names(a)
        a$gdppc <- exp(a$lgdppc)
        a$bpc <- exp(a$lbpc)
        a$present08 <- a$present08/1000
        a$Area <- a$Area/1000
        
        
        cont_test <- c("outcome","seats08","t.mayor",
                       "qe","gdppc","bpc","present08","poverty","Area","rent.inc","pt_base") 
        
        a %>% select(all_of(cont_test))
        
        
        jx <- c(3,3,3,3,3,3,3,3,3,3,3,3,3)
        
        a1 <- subset(a, MainSample == 1)
        a0 <- filter(a, MainSample == 0 & !tse %in% a1$tse)
        a0 <- unique(a0, by="tse")
        
        res <- data.frame(matrix(0,length(cont_test)+1,3))
        
        i <- 1
        
        for (i in 1:length(cont_test)){
          
          vec1 <- pull(a1, cont_test[i]) 
          vec0 <-  pull(a0, cont_test[i]) 
          pv <- t.test(vec1,vec0)$p.value
          
          stars <- ""
          if(pv <= 0.1) stars <- "+"
          if(pv <= 0.05) stars <- "*"
          
          res[i,1] <- format(round(mean(vec1),jx[i]),nsmall=jx[i])
          res[i,2] <- format(round(mean(vec0),jx[i]),nsmall=jx[i])
          res[i,3] <- paste0(format(round(mean(vec1)-mean(vec0),jx[i]),nsmall=jx[i]),stars)
          
          
          
        } 
        
        
        res[nrow(res),1] <- nrow(a1)
        res[nrow(res),2] <- nrow(a0)
        res[nrow(res),3] <-nrow(a1)-nrow(a0)
        
        nami <- c("Council Expansion","Seats in 2008","Mayor's Support",
                  "Votes per Seat","GDP pc","Budget pc","Voters", "Poverty","Area",
                  "Mayor's party in PMDB/PP/PL/PTB",
                  "Mayor's party is NOT PSDB/PPS/DEM")
        
        spaces <- nrow(res)
        nam <- c(nami,"Municipalities")
        tt <- .LatexTable(res,nam,spaces)
        print(tt, row.names=F)
  
        
    #===================================================================================================================
        
        
    # Table A.7 =========================================================================================================
        
        
        setwd("C:/Users/ander/Dropbox/Papers/Legsize/replication_note/")
        m <- fread("data_reelection.csv")
        m$year <- ifelse(m$election==2012,1,0)
        mx <- filter(m, run==1)
        
        
        tab <- data.frame(matrix(0,15,6))
        for (i in 1:6){
          
          if(i==1) reg <- felm(run ~ (year + change * year)  | 0 | 0 | tse, data=m )
          if(i==2) reg <- felm(run ~ (year + change : year)  | tse | 0 | tse, data=m )
          if(i==3) reg <- felm(run ~ lead*(year + change : year) | tse | 0 | tse, data=m )
          if(i==4) reg <- felm(win ~ (year + change * year) | 0 | 0 | tse, data=mx)
          if(i==5) reg <- felm(win ~ (year + change : year)  | tse | 0 | tse, data=mx )
          if(i==6) reg <- felm(win ~ lead*(year + change : year)  | tse | 0 | tse, data=mx )
          
          r <- data.frame(summary(reg)$coefficients)
          names(r)[3] <- "tt"
          r$star <- "*"
          r$star[abs(r$tt)<1.96] <- "+"
          r$star[abs(r$tt)<1.645] <- " "
          
          if(i %in% c(1,4)) {tab[1,i] <- paste0(round(r["(Intercept)",1],3),r["(Intercept)",5])}else{tab[1,i] <- ""}
          if(i %in% c(1,4)) {tab[2,i] <- paste0("(",format(round(r["(Intercept)",2],3),nsmall=3),")")}else{tab[2,i] <- ""}
          
          if(i %in% c(1,4)) {tab[3,i] <- paste0(round(r["change",1],3),r["change",5])}else{tab[3,i] <- " "}
          if(i %in% c(1,4)) {tab[4,i] <- paste0("(",format(round(r["change",2],3),nsmall=3),")")}else{tab[4,i] <- " "}
          
          
          tab[5,i] <- paste0(format(round(r["year",1],3),nsmall=3),r["year",5])
          tab[6,i] <- paste0("(",format(round(r["year",2],3),nsmall=3),")")
          tab[7,i] <-  paste0(format(round(r["year:change",1],3),nsmall=3),r["year:change",5])
          tab[8,i] <- paste0("(",format(round(r["year:change",2],3),nsmall=3),")")
          
          if(i %in% c(3,6)){
            
            tab[9,i] <-  paste0(format(round(r["lead",1],3),nsmall=3),r["lead",5])
            tab[10,i] <- paste0("(",format(round(r["lead",2],3),nsmall=3),")")
            
            tab[11,i] <-  paste0(format(round(r["lead:year",1],3),nsmall=3),r["lead:year",5])
            tab[12,i] <- paste0("(",format(round(r["lead:year",2],3),nsmall=3),")")
            
            tab[13,i] <-  paste0(format(round(r["lead:year:change",1],3),nsmall=3),r["lead:year:change",5])
            tab[14,i] <- paste0("(",format(round(r["lead:year:change",2],3),nsmall=3),")")
          }
          
          if(i %in% c(1:3)) {tab[15,i] <- nrow(m)}else{tab[15,i] <- nrow(mx)} 
        }
        
        tab[is.na(tab) | tab==0] <- " "
        tab
        
        
        # LATEX TABLE
        
        # names of columns
        lista_rownames <- c("Intercept"," ",
                            "Seat Increase (A)"," ",
                            "2012 Election (B)"," ",
                            "A x B"," ",
                            "Coalition (C)"," ",
                            "A x B "," ",
                            "A x B x C"," ",
                            "Observations")
        # spaces after the following lines
        spaces <- c(12)
        
        table <- .LatexTable(tab,lista_rownames,spaces)
        print(table, row.names=FALSE)
        
    #===================================================================================================================
        
        
    # Figure A.1 =======================================================================================================
        
        
        new <- filter(base, MainSample==1)  
        band <- 1.7
        data <- data.frame(subset(new, abs(run) <= band   )  )
        szrd <- 12
        ll <- -.4
        lh <-  .4
        bins <- 14
        
        data$nvar <- data$outcome
        data$var <- summary(felm( nvar~1 |fe|0|0 ,data=data) )$res
        brr <- round(seq(from=-band, to=band, length=7),2)
        prdd1 <- .RDPlot(data=data,cut=0,band=band,bins=bins,degree=1,
                         x.title="Running Variable",
                         y.title="Normalized outcome",
                         title="DV: Council Expansion (1=YES)",
                         con=0.95,sz=szrd,brk=brr,
                         limi1=ll,limi2=lh) 
    
        
        data$nvar <- data$t.mayor
        data$var <- summary(felm( nvar~1 |fe|0|0 ,data=data) )$res
        
        brr <- round(seq(from=-band, to=band, length=7),2)
        prdd2 <- .RDPlot(data=data,cut=0,band=band,bins=bins,degree=1,
                         x.title="Running Variable",
                         y.title="Normalized outcome",
                         title="DV: Mayor's Support",
                         con=0.95,sz=szrd,brk=brr,
                         limi1=-.2,limi2=.2) 
  
        
        grid.arrange(prdd1, prdd2,   ncol=2)
        rm(new,data,prdd2,prdd1,brr,band,szrd,ll,lh,bins)
        
        
    #===================================================================================================================
        
        
    # Figure A.2 =======================================================================================================
        
        # keeps only observations where the close election was: coalition vs. opposition
        new <- filter(base, MainSample==1)    
      
        for (out in 1:2){
          
          if(out==1){hlimi <- 0.25}else{hlimi <- 0.4}
          if(out==1){llimi <- -1}else{llimi <- 0}
          
          if(out==1){new$var <- new$outcome}else{new$var <- new$t.mayor}
          bandvector <- seq(from=0.75, to=1.5, by=0.25)
          res <- data.table(matrix(0,length(bandvector),7))
          
          ddd <- new %>% select(fe)
          datafe <- dummy_cols(ddd)
          datafe <- datafe[,3:ncol(datafe)]
          datacon <- new %>% select(all_of(controls))
          ddd <- new %>% select(reg)
          datasta <- dummy_cols(ddd)
          datasta <- datasta[,3:ncol(datasta)]
          
          b <- rdbwselect(new$var, new$run, p=1, q=qq, bwselect = btype)$bws ; b
          
          
          for (i in 1:3){
            
            if(i == 1) cv <- datafe
            if(i == 2) cv <- cbind(datafe,datacon)
            if(i == 3) cv <- cbind(datafe,datacon,datasta)
            
            for (j in 1:length(bandvector)){
              
              b1 <-  b[1]*bandvector[j] 
              b2 <-  b1*(b[3]/b[1])
              
              reg <- rdrobust(new$var,new$run, covs=cv, b=b2, h=b1, q=qq, vce=vart) 
              reg2 <- rdrobust(new$var,new$run, covs=cv, b=b2, h=b1, q=qq, level=90, vce=vart)
              #summary(reg)
              typeci <- 3
              
              res[j,1] <- reg$coef[1]
              res[j,2] <- bandvector[j]*100
              res[j,3] <- reg$ci[typeci,1]
              res[j,4] <- reg$ci[typeci,2]
              res[j,5] <- reg2$ci[typeci,1]
              res[j,6] <- reg2$ci[typeci,2]
              res[j,7] <- b1
            }
            
            print(i)
            if(out==1 & i == 1) p.title <- "COUNCIL EXPANSION"
            if(out==1 & i == 2) p.title <- "COUNCIL EXPANSION (Covariates)"
            if(out==1 & i == 3) p.title <- "COUNCIL EXPANSION (Covariates and State F.E.)"
            if(out==2 & i == 1) p.title <- "MAYOR'S SUPPORT"
            if(out==2 & i == 2) p.title <- "MAYOR'S SUPPORT (Covariates)"
            if(out==2 & i == 3) p.title <- "MAYOR'S SUPPORT (Covariates and State F.E.)"
            
            
            names(res) <- c("co","band","ci95a","ci95b","ci90a","ci90b","band2")
            pl <- ggplot(data=res, aes(x=band, y=co))   +    theme_bw() +
              ylim(llimi,hlimi)+
              theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
              geom_hline(yintercept = 0 , lty=3, alpha=.75) +
              theme(plot.title = element_text(family = "A", size=13))+
              xlab("% of Optimal Bandwidth") + ylab("RDD Effect") +
              theme(axis.text = element_text(colour="black", size=13, family="A")) +
              #theme(axis.title = element_text(colour="black", size=13, family="A")) +
              theme(axis.title = element_blank()) +
              theme(plot.title = element_text(colour="black", size=11, family="A")) +
              geom_errorbar(aes(x=band, ymin=ci95a, ymax=ci95b), colour="gray15",  linewidth=1, alpha=.75, width=0) +
              geom_errorbar(aes(x=band, ymin=ci90a, ymax=ci90b), colour="gray25",  linewidth=3, alpha=.5, width=0) +
              geom_point(aes(x=band, y=co), colour="gray15",  size=4, alpha=.75) +
              scale_x_continuous(breaks=bandvector*100)+
              ggtitle(p.title)
            
            if(out==1 & i == 1) p11 <- pl
            if(out==1 & i == 2) p12 <- pl
            if(out==1 & i == 3) p13 <- pl
            if(out==2 & i == 1) p21 <- pl
            if(out==2 & i == 2) p22 <- pl
            if(out==2 & i == 3) p23 <- pl
            
            
          }
          
          
        }
        
        grid.arrange(p11,p21,p12,p22,p13,p23,nrow=3, 
                     bottom="% of the Optimal Bandwith", left="RDD Effect")
        
    #===================================================================================================================
        
        
    # Figure A.3 =======================================================================================================
        
       
        controls <- c("qe","lgdppc","lbpc","lvoters","lpoverty","rent.inc","lArea","pt_base")  #lbudget
        
        # keeps only observations where the close election was: coalition vs. opposition
        new <- data.table(filter(base, MainSample==1)) 
     
        for (out in 1:2){
          
          if(out==1){new$var <- new$outcome}else{new$var <- new$t.mayor}
          res <- data.frame(matrix(0,5,3))
          bandvector <- seq(from=0.75, to=1.5, by=.25)
          res <- data.table(matrix(0,length(bandvector),7))
          if(out==1){hlimi <- 0.25}else{hlimi <- 0.4}
          if(out==1){llimi <- -1}else{llimi <- 0}
          
          ddd <- new %>% select(fe)
          datafe <- dummy_cols(ddd)
          datafe <- datafe[,3:ncol(datafe)]
          datacon <- new %>% select(all_of(controls))
          ddd <- new %>% select(reg)
          datasta <- dummy_cols(ddd)
          datasta <- datasta[,3:ncol(datasta)]
          
          b <- rdbwselect(new$var, new$run, c = 0, p=2, q=5, bwselect = btype)$bws ; b/sdd
          
          
          for (i in 1:3){
            
            if(i == 1) cv <- datafe
            if(i == 2) cv <- cbind(datafe,datacon)
            if(i == 3) cv <- cbind(datafe,datacon,datasta)
            
            
            for (j in 1:length(bandvector)){
              
              b1 <- b[1]*bandvector[j]
              b2 <- b1*(b[3]/b[1])
              
              reg <-  rdrobust(new$var,new$run, covs=cv, b=b2, h=b1,  q=5,p =2,           vce=vart) 
              reg2 <- rdrobust(new$var,new$run, covs=cv, b=b2, h=b1,  q=5,p =2, level=90, vce=vart)
              
              res[j,1] <- reg$coef[1]
              res[j,2] <- bandvector[j]*100
              res[j,3] <- reg$ci[3,1]
              res[j,4] <- reg$ci[3,2]
              res[j,5] <- reg2$ci[3,1]
              res[j,6] <- reg2$ci[3,2]
              res[j,7] <- b1/sdd
            }
            
            print(i)
            if(out==1 & i == 1) p.title <- "COUNCIL EXPANSION"
            if(out==1 & i == 2) p.title <- "COUNCIL EXPANSION (Covariates)"
            if(out==1 & i == 3) p.title <- "COUNCIL EXPANSION (Covariates and State F.E.)"
            if(out==2 & i == 1) p.title <- "MAYOR'S SUPPORT"
            if(out==2 & i == 2) p.title <- "MAYOR'S SUPPORT (Covariates)"
            if(out==2 & i == 3) p.title <- "MAYOR'S SUPPORT (Covariates and State F.E.)"
            
            
            names(res) <- c("co","band","ci95a","ci95b","ci90a","ci90b","b")
            pl <- ggplot(data=res, aes(x=band, y=co))   +    theme_bw() +
              ylim(llimi,hlimi)+
              theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
              geom_hline(yintercept = 0 , lty=3, alpha=.75) +
              theme(plot.title = element_text(family = "A", size=13))+
              xlab("% of Optimal Bandwidth") + ylab("RDD Effect") +
              theme(axis.text = element_text(colour="black", size=13, family="A")) +
              #theme(axis.title = element_text(colour="black", size=13, family="A")) +
              theme(axis.title = element_blank()) +
              theme(plot.title = element_text(colour="black", size=11, family="A")) +
              geom_errorbar(aes(x=band, ymin=ci95a, ymax=ci95b), colour="gray15",  linewidth=1, alpha=.75, width=0) +
              geom_errorbar(aes(x=band, ymin=ci90a, ymax=ci90b), colour="gray25",  linewidth=3, alpha=.5, width=0) +
              geom_point(aes(x=band, y=co), colour="gray15",  size=4, alpha=.75) +
              scale_x_continuous(breaks=bandvector*100)+
              ggtitle(p.title)
            
            if(out==1 & i == 1) p11 <- pl
            if(out==1 & i == 2) p12 <- pl
            if(out==1 & i == 3) p13 <- pl
            if(out==2 & i == 1) p21 <- pl
            if(out==2 & i == 2) p22 <- pl
            if(out==2 & i == 3) p23 <- pl
            
            
          }
          
          
        }
        
        grid.arrange(p11,p21,p12,p22,p13,p23,nrow=3, bottom="% of the Optimal Bandwith", left="RDD Effect")
        
        
        
    #===================================================================================================================
        
        
    # Figure A.4 =======================================================================================================
        
        # keeps only observations where the close election was: coalition vs. opposition
        new <- data.table( filter(base, MainSample==1 & !is.na(share))  )
        new$var <- new$outcome
        
        # Keep only data with donor information
        new$patr <- ifelse(new$share2 > median(new$share2),1,0)
        
        ddd <- new %>% select(fe)
        datafex <- dummy_cols(ddd)
        datafex <- datafex[,3:ncol(datafex)]
        dataconx <- new %>% select(all_of(controls))
        ddd <- new %>% select(reg)
        dataregx <- dummy_cols(ddd)
        dataregx <- dataregx[,3:ncol(dataregx)]
     
        
        # Optimal bandwidth
        b <- rdbwselect(new$var, new$run, p=1, q=qq, bwselect = btype)$bws 
        b1 <- b[1]
        b2 <- b1*(b[3]/b[1])
        
        res <- data.frame(matrix(0,3,7))
        
        bb <- new
        
        # Calculate individual regressions for each subsample based on patronage
        for (i in 0:1){
          
          bb$i <- abs(i-bb$patr)
          bb$tr <- bb$t*bb$run
          bb$ti <- bb$t*bb$i
          bb$tir <- bb$t*bb$i*bb$run
          bb$ir <- bb$i*bb$run
          cv <- cbind(bb$i,bb$ir,bb$tir,bb$ti,datafex,dataconx,dataregx)
          
          reg <- rdrobust(bb$var,bb$run,  b=b2, h=b1,covs=cv  , masspoints="off", vce=vart) 
          reg2 <- rdrobust(bb$var,bb$run,  b=b2, h=b1,covs=cv , level=90 , masspoints="off", vce=vart) 
 
          stars <- ""
          if(reg2$ci[3,2]*reg2$ci[3,1]>0) stars <- "+"
          if(reg$ci[3,2]*reg$ci[3,1]>0) stars <- "*"
          
          
          res[1+i,1] <- reg$coef[1]
          res[1+i,2] <- reg$ci[3,1]
          res[1+i,3] <- reg$ci[3,2]
          res[1+i,4] <- reg2$ci[3,1]
          res[1+i,5] <- reg2$ci[3,2]
          res[1+i,6] <- stars
          res[1+i,7] <- b1
          
        }
        
        res[3,1] <- res[2,1] - res[1,1]
        
        # Bootstrap the CI for the difference in coefficients
        set.seed(14618)
        results <- boot(data=new,statistic=.DiffRDD , R=500,
                        b1=b1,b2=b2, 
                        controls=controls,variable="patr",
                        parallel = "snow", ncpus = 4)
        
        ci = boot.ci(results, type=c("perc"), conf=.95) ;  ci 
        ci2 = boot.ci(results, type=c("perc"), conf=.90) ;  ci2
        res[3,2] <-    ci$perc[4]
        res[3,3] <-    ci$perc[5]
        stars <- ""
        if(ci2$perc[4]*ci2$perc[5]>0) stars <- "+"
        if(ci$perc[4]*ci$perc[5]>0) stars <- "*"
        res[3,6] <- stars
        names(res) <- c("co","lo","hi","lo2","hi2","stars","band")
        res$nm <- c("Low \n Patronage","High \n Patronage","Difference")
        res$x <- as.factor(1:3)
        res_patr <- res
        
        ph1x <- ggplot(data=res_patr, aes(x=x, y=co))   + theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_hline(yintercept=0, lty=3)+
          geom_point( size=5)+
          geom_errorbar(aes(ymin=lo,ymax=hi),width=0.0, linewidth=1.25)+
          scale_x_discrete(labels=res_patr$nm) +
          theme(axis.text = element_text(colour="black", size=13, family="A")) +
          theme(axis.title = element_blank())+
          ylab("Coefficient (RD Effect)") 
          
        
        # Restrict the data to observations within the optimal band
        data_short <- filter(new, abs(run) < b1)
        data_short$var <- data_short$outcome
        data_short$wght <- b1 - abs(data_short$run)
        
        formula2 <- as.formula( paste0( "var ~ t*run*share2+", paste(controls,collapse="+") ,"|reg+fe|0|0"  ))
        reg <- felm( formula2,  weights=data_short$wght, data=data_short) ; summary(reg, robust=TRUE)
        
        rr <- summary(reg)$coefficients
        vv <- vcov(reg)
        
        vec <- round(seq(from= -2.5, to = 1.5,length.out=5),2)
        vec2 <- round(seq(from= -1, to = 1,length.out=5),1) 
        ph2x <- .NewPlot(data=data_short,r=rr,v=vv,xvar="share2",minD=-1,maxD=1,bins1=12,bins2=14,
                         shift=2.5,bb=vec,degree=1,multi=12,oneside=0,
                         x.title="Patronage Level",y.title="",p.title="",bb3=vec2,pt=0) 
        
        
        grid.arrange(ph1x,ph2x,ncol=2,widths=c(0.4, 0.6),
                     left = text_grob("RDD Effect", size = 13,  family="A", rot = 90)) 
        
        rm(bb,new,rr,vv,data_short,res_patr,res,dataregx,dataconx,datafex)
        
        
    #===================================================================================================================
      
        
    # Figure A.5 =======================================================================================================
        
        
        # keeps only observations where the close election was: coalition vs. opposition
        new <- filter(base, MainSample==1 )    
        new$var <- new$outcome
        new$lp <- as.numeric(substring(new$no_win,1,2))*new$last.w + 
          as.numeric(substring(new$op.no_win,1,2))*new$last.l
        new$rent.council <- ifelse(new$lp %in% c(15,11,22,14),1,0)
       
        # Optimal bandwidth
        b <- rdbwselect(new$var, new$run, p=1, q=qq, bwselect = btype)$bws ; b
        b1 <- b[1]
        b2 <- b1*(b[3]/b[1])
        
        
        ddd <- new %>% select(fe)
        datafex <- dummy_cols(ddd)
        datafex <- datafex[,3:ncol(datafex)]
        dataconx <- new %>% select(all_of(controls))
        ddd <- new %>% select(reg)
        dataregx <- dummy_cols(ddd)
        dataregx <- dataregx[,3:ncol(dataregx)]
        
        
        res <- data.frame(matrix(0,3,7))
        
        bb <- new
        
        # Calculate individual regressions for each subsample based on patronage
        for (i in 0:1){
          
          bb$i <- abs(i-bb$rent.council)
          bb$tr <- bb$t*bb$run
          bb$ti <- bb$t*bb$i
          bb$tir <- bb$t*bb$i*bb$run
          bb$ir <- bb$i*bb$run
          cv <- cbind(bb$i,bb$ir,bb$tir,bb$ti,datafex,dataconx,dataregx)
          
          reg <- rdrobust(bb$var,bb$run,  b=b2,  h=b1,covs=cv ,q=qq, masspoints="off", vce=vart) 
          reg2 <- rdrobust(bb$var,bb$run,  b=b2, h=b1,covs=cv ,q=qq,  level=90 , masspoints="off", vce=vart) 
           
          stars <- ""
          if(reg2$ci[3,2]*reg2$ci[3,1]>0) stars <- "+"
          if(reg$ci[3,2]*reg$ci[3,1]>0) stars <- "*"
          
          
          res[1+i,1] <- reg$coef[1]
          res[1+i,2] <- reg$ci[3,1]
          res[1+i,3] <- reg$ci[3,2]
          res[1+i,4] <- reg2$ci[3,1]
          res[1+i,5] <- reg2$ci[3,2]
          res[1+i,6] <- stars
          res[1+i,7] <- b1
          
        }
        
        res[3,1] <- res[2,1] - res[1,1]
        
        # Bootstrap the CI for the difference in coefficients
        set.seed(14618)
        results <- boot(data=new,statistic=.DiffRDD , R=500,
                        b1=b1,b2=b2, 
                        controls=controls,variable="rent.council",
                        parallel = "snow", ncpus = 6)
        
        ci = boot.ci(results, type=c("perc","basic"), conf=.95) ;  ci 
        ci2 = boot.ci(results, type=c("perc","basic"), conf=.90) ;  ci2
        res[3,2] <-    ci$basic[4]
        res[3,3] <-    ci$basic[5]
        stars <- ""
        if(ci2$basic[4]*ci2$basic[5]>0) stars <- "+"
        if(ci$basic[4]*ci$basic[5]>0) stars <- "*"
        res[3,6] <- stars
        names(res) <- c("co","lo","hi","lo2","hi2","stars","band")
        res$nm <- c("Other \n Parties","PMDB,PTB, \n PP,PL","Difference")
        res$x <- as.factor(1:3)
        res_part <- res
        
        ph5 <- ggplot(data=res_part, aes(x=x, y=co))   + theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_hline(yintercept=0, lty=3)+
          geom_point( size=5)+
          geom_errorbar(aes(ymin=lo,ymax=hi),width=0.0, linewidth=1.25)+
          scale_x_discrete(labels=res_part$nm) +
          theme(axis.text = element_text(colour="black", size=13, family="A")) +
          theme(plot.title = element_text(colour="black", size=13, family="A")) +
          theme(axis.title = element_blank())+  ggtitle("IV: Councilor's Partisanship") +
          ylab("Coefficient (RD Effect)") 

        
        # Restrict the data to observations within the optimal band
        data_short <- filter(new, abs(run) < b1)
        data_short$var <- data_short$outcome
        data_short$wght <- b1 - abs(data_short$run)
        data_short$v <- data_short$v/100
 
        formula2 <- as.formula( paste0( "var ~ t*run*v+", paste(controls,collapse="+") ,"|reg+fe|0|ibge"  ))
        reg <- felm( formula2,  weights=data_short$wght, data=data_short) ; summary(reg)
        
        
        rr <- summary(reg)$coefficients
        vv <- vcov(reg)
        
        vec <- round(seq(from= -3, to = 2,length.out=6),2)
        vec2 <- round(seq(from= 0, to = .3,length.out=5),1) 
        ph4 <- .NewPlot(data=data_short,r=rr,v=vv,xvar="v",minD=0,maxD=.3,bins1=12,bins2=14,
                        shift=3,bb=vec,degree=1,multi=6,oneside=0,
                        x.title="Last Councilor's Vote as Share of Coalition",y.title="",
                        p.title="IV: Councilor's Safety",bb3=vec2,pt=1) 
        
        grid.arrange(ph5,ph4,ncol=2,widths=c(0.5, 0.5),left = text_grob("RDD Effect", size = 13,  family="A", rot = 90)) 
  
        
    #===================================================================================================================
       
        
    # Figure A.6 ========================================================================================================
        
        new <- filter(base, MainSample==1)
        test1 <- rddensity(X=new$run,c=0,p=1,q=6)
        summary(test1)
        rdplotdensity(test1, X=new$run, histFillCol="gray50",CIcol=c("gray75","gray25"),
                      lcol=c("gray75","gray25"), plotRange=c(-3,3))
       

    #===================================================================================================================
        
        
        
        
        
  
        